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Abstract 



Ensemble Kalman filter techniques are widely used to assimilate observations into dy- 
£\j ■ namical models. The phase space dimension is typically much larger than the number of 

ensemble members which leads to inaccurate results in the computed covariance matri- 
| ces. These inaccuracies can lead, among other things, to spurious long range correlations 

which can be eliminated by Schur-product-based localization techniques. In this paper, 
we propose a new technique for implementing such localization techniques within the class 
of ensemble transform/square root Kalman filters. Our approach relies on a continuous 
embedding of the Kalman filter update for the ensemble members, i.e., we state an ordi- 
nary differential equation (ODE) whose solutions, over a unit time interval, are equivalent 
to the Kalman filter update. The ODE formulation forms a gradient system with the ob- 
servations as a cost functional. Besides localization, the new ODE ensemble formulation 
qq ■ should also find useful applications in the context of nonlinear observation operators and 

C — ■ . observations arriving continuously in time. 

Keywords. Data assimilation, ensemble Kalman filter, localization, continuous Kalman filter 

G\ ■ 

o : 

O ■ 1 Introduction 

> : 

We consider ordinary differential equations 

.5?: *=/(«,*) a) 

with state variable xel" Initial conditions at time to are not precisely known and we assume 
instead that 

x(t ) ~N(xo,B), (2) 

where N(x , B) denotes an n-dimensional Gaussian distribution with mean x G M n and co- 
variance matrix B G M nxn . We also assume that we obtain measurements y(ii) G R fc at 
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discrete times tj > to, j — 0,1, ... , M, subject to measurement errors, which are also Gaussian 
distributed with zero mean and covariance matrix R G 



ikxk 



i.e. 



y(t,)-Hx(t,)~N(0,R). 



(3) 



Here H G IR fcxn is the (linear) measurement operator. 

Data assimilation is the task to combine the model ([T]) (here assumed to be perfect), the 
knowledge about the initial conditions (j2J) and available measurements (ES) in a predic t ion of the 
probability distribution of the solution at any time t > to. We refer to iLewis et al.l (120061 ) for 
a detailed introduction and available approaches to data assimilation. In this paper, we focus 
on the ensemble Kalman filter (EnKF) method, originally proposed by Evensen fsee lEyensen 
(120061 ) for a recent ac count) and, in p articular, on ensemble transform ( Bishop et all 12001 ). 

), and ensemble square root filters ( Tippct t et al.l . 120031 ) 



ensemble adjustment ( Anderson] . 200 



and their sequential implementation (I Whitaker and Hamilll . 120021 ; lAndersonl . 120031 ) 



The EnKF relies on the simultaneous propagation of m independent solutions Xj(t), i 
1, . . . , m, from which we can extract an empirical mean 



^ m 

x(t) = -$>,(£) 
i=l 



(4) 



and an empirical covariance matrix 

1 m 

p(t) = — - £ (Mt) - x(t)) (Mt) - x(t)) T . 



(5) 



1=1 



In typical applications from meteorology, the ensemble size m is much smaller than the dimen- 
sion n of the model phase space and, more importantly, also much smaller than the number 
of positive Lyapunov exponents. Hence P(t) is highly rank deficie nt which can lead to un- 
reliab le pre dictions. Ensemble localization has been introduced by iHoutekamer and Mitchell 
(120011 ) and lHamill et al.l ( 1200 ll ) to overcome this problem. However, only two techniques are 
currently available to implement Schur-product-based localization within the framework of en- 
semble transform/square r oot Kalman filters. The first option i s prov ided by a sequential 
processing of observations (jWhitaker and Hamilll. |2002|; lAndersonl. 120031 ). while the determin- 
istic ensemble Kalman filter (DEnKF) of ISakov and Okd (j2008af ) is a second, more recent, 
option. The DEnKF results in an approximate implementation of e nsemble t ransform/square 



root Kalman filters. W e also mention box/local analysis methods (lEvensenl . 120031 ; lOtt et al. 



2004 ; iHunt et al.l . 120071 1 . which assimilate data locally in physical space and which therefore 



possess a "built in" localization. 

In this paper, we demonstrate that techniques proposed by Bergemann et all ( 2009 ) for the 
filter analysis step can be further generalized to an ordinary differential equation (ODE) formu- 
lation in terms of the ensemble members Xj, i — 1, . . . , m. This formulation is subsequently used 
to derive an easy to implement localized ensemble Kalman filter, which can process observations 
simultaneously and can be extended to nonlinear observation operators. 
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2 Background material 



We su mmarize a numb er of key results and techniques regarding ensemble Kalman filters. We 
refer to lEvensenl (120061 ) for an introduction and in-depth discussion of such filters. 



2.1 Kalman analysis step 

Let n denote the dimension of the phase space of the problem. We consider an ensemble of m 
members Xj(i) G W 1 which we collect in a matrix X(t) G M nXTn . In terms of X, the ensemble 
mean is given by 

(6) 



x(t) = — X(t)e G R 
m 

and we introduce the ensemble deviation matrix 

X'(t) = X(t) -x(£)e T G 

1,...,1) T G M m . 



(7) 



where e = (1, . 

We now describe the basic Kalman analysis step. Let x/- and X^- denote the forecast mean 
and deviation matrix, respectively. The ensemble mean is updated according to 



K (Hx 



where 



K = P/H T (HP/H T + R) 1 

is the Kalman gain matrix with empirical covariance matrix 

1 



m — 1 



-X' X 



(8) 
(9) 

(10) 



and R G R kxk is the measurement error covariance matrix. 

While the update of the mean is common to most ensemble Kalman filters, the update of 
the ensemble deviation matrix X^ can be implemented in several ways. In this paper, we focus 
on ensemble update techniques that employ either a transformation of the form 



X' = AX' 



/ 



with an appropriate matrix A G 



(lAndersonl . l200ll ) or a transformation 
K = X/T 



(11) 



(12) 



with an appropriate T G M. mxm flBishop et al!l200ll : IWhitaker and Hamuli 120021 : iTippett et all 
20031 : lEvensenl . 120041 ) . The matrices A and T are chosen such that the resulting ensemble 
deviation matrix X' satisfies 



m — 1 



-X'fXl 



(I-KH)^. 



(13) 
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It has been shown by iTippett et al.l (120031 ) that both formulations ffTTl) and ( I12p can be made 
equivalent not only in terms of (fT3"j) but also in terms of the resulting ensemble deviation matrix 
X' a . Since n > m in most applica tions, formulatio n ( lT2"j) is generally preferred except when 
working in a sequential framework (lAndersonl . 120031 1 . 

Note that the transformation matrix T should also satisfy Te = e to guarantee X^e = 
( jWang et all . I2004J ; Livings et afl 120081 : ISakov and Okel . l2008bl ). Otherwise, the update of the 
ensemble deviation matrix would affect the update of the ensemble mean. 



S everal methods h ave been proposed 
fl2009h ) that satisfy 



and IBer geman n et al , _ 

Sakov and Okel ( 2008a ) suggest to use ffTTl) with 



recently (including ISakov and Okel (j2008al ) 
T3l) only approximately. More specifically, 



I 



^KH, 

2 



(14) 



while IBergemann et al.l ( 120091 ) use numerical approximations to the underlying ODE formula- 
tion 

d 1 



ds 



Y 



2m -2 



YY T H T R X HY 



(15) 



in a fictitious time s G [0, 1] See, for example, ISimonl ( 120061 ) for a derivation of (jl5|) . The 
initial condition is Y(0) = and the updated ensemble deviation matrix, which satisfies (1131) 
exactly, is provided by the solution at time s — 1, i.e. 

X' = Y(1). (16) 

A typica l numerical implementat ion of ( ITS"]) uses two or four time-steps with the forward Euler 
method ( Bergemann et aP . 2009 ). The resulting transformation of the forecast into the ana- 
lyzed ensemble deviation matrix is of the form ( TTTj) with A defined through the time-stepping 
method. 

Note that the Kalman gain matrix (0 is equivalent to 

K = P a H T R-\ (17) 

which is advantageous in connection with (fT5|) since only the inversion of the measurement 
error covariance matrix R G R kxk is now required to implement a complete Kalman analysis 
step. Algorithmically, one would first update the ensemble deviation matrix using ( TT5"j) . then 
form the analysed ensemble covariance matrix P a = X^[X^] T /(m — 1) as well as the Kalman 
gain matrix (fTTf) . and finally update the ensemble mean using (JH]). 

All methods discussed so far have in common that the Kalman update increments for the 
ensemble mean and ensemble deviation matrix lie in a m — 1 dimensional subspace, denoted by 
§j C R n . This space is defin ed by the range/image of the forecast ensemble deviation matrix 



Xj. IBergemann et al.l (120091 ) introduced a continuous matrix factorization algorithm for the 



ensemble X(t), which automatically produces orthog onal vectors that span §f . 

It is common practice to apply variance inflation (lAnderson and Anderson! . 119991 ) to X'(^) 
before the forecasted ensemble is updated under the Kalman filter analysis step, i.e., the Kalman 
analysis step uses 



X 



x{tj)e T + 6X'{t j ), 



(18) 



where 5 > 1 is an inflation factor, instead of Xj = X(tj) 
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2.2 Localization 



The idea of localization, as proposed by iHoutekamer and Mitchell! ( 120011 ). is to replace the 
matrices HP f and HP /H T in the Kalman gain matrix (|9]) by 



BF f = C loc>1 o (HP f ) , HP f H T = C loCi2 o (HP / H T ) , 



(19) 



pnxk 



and C 



lnr„? 



okxk 



are a ppropriate localization matrices based 



respectively, where Ci oc ,i G 
on filter functions suggested by (iGaspari and Cohnl . Il999l ) and CoY denotes the Schur product 
of two matrices C and Y of identical dimension, i.e., 

(C O Y)y = (C)y (Y)y (20) 

for all indices We denote the resulting modified Kalman gain matrix by Ki oc f, i.e., 

K Ioc>f = (HP f ) T (hP^IF + R) ' . (21) 



Localization was also proposed by lHamill et al. (l200lh with the only difference that localization 
is not applied to HP/H T . 

Alternatively, one can localize the Kalman gain matrix formulation (fTTj) and use 



K 



loc,a 



fHP„) T R 



(22) 



instead of in an ensemble Kalman filter. Note that (1211) and (1221) are not equivalent in 
general and that formulation ([22]) is easier to implement. 

Based on these modified Kalman gain matrices, a Schur-product-based localization is easy 
to apply to the update (El) of the ensemble m ean and to ensemble deviation updates that use 
pe rturbed observations ([Burgers et all 119981). which i s esse ntially the localization approach 



of IHoutekamer and Mitchell! (120011 ) and lHamill et al.l (120011 ). However, the popular class of 
ensemble transform/square root filters, based on (TT2|) . has not yet been amenable to Schur- 
product-based localizations e xcept when observations are treated sequentially, i.e., when k = 1 
in each transformation step ( Whitaker and Hamilll . 2002 ). 

It is feasible that localizations can be implemented for ensemble deviation updates of 
the form (1111) th rough an appropriate modification of the ensemble adjustment technique of 
Anderson! ( 200ll ). However, such a modification would lead to a computationally expensive 
implementation of Sch ur-product-based localizations. The recently proposed DEnKF filter of 
Sakov and Oke (j2008al ). on the other hand, leads to a computationally feasible implementation 
with the localization directly applied to (THj) . i.e., one uses 



A 



I — -Ki OCj fH 



(23) 



in Am 

We note that localization implies in general that the update increments for the ensem- 
ble mean and the ensemble deviation matrix no longer lie in the subspace §y defined by the 
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range/image of X^. While this is a desirable property on the one hand, it can lead to unbal- 
anced fields in the analyzed ensemble X„ on the other hand. T his has been investigated, for 
example, by iHoutekamer and Mitchelll (120051 ) and iKepertl (120091 ). 

We finally mention an a lterna t ive approach to lo calization. The box/local EnKF filters of 



Evensenl (120031 ); lOtt et al.l (120041 ); iHunt et al.l (120071 ) assimilate data locally in physical space 
and possess a "built in" localization based on the spatial structure of the underlying partial 
differential equation model. 



3 Localization based on continuous ensemble updates 

We now describe an alternative for introducing localization, which is based on a generalization 
of the ODE formulation (I15p and which leads to an ODE formulation directly in the ensemble 
members Xj. 

We first note that the Kalman update ([8]) for the ensemble mean can also be formulated in 
terms of an ODE, i.e., 

—x = — YY T H T R _1 (Hx - y) (24) 

ds m - 1 v ' y 1 

with x(0) = x/ and x a = x(l). See, for example, ISimonl (120061 ) for a derivation of ( 12~4"|) . 

To further reveal the underlying mathematical structure of ( 1151) and (1241) . we introduce the 

cost functional 

5(x) = l - (Hx - y) r R 1 (Hx - y) (25) 

for each set of observations. Next we combine ( 1151) and (12~4"|) to give rise to the differential 
equations 

^x, = -ip {V Xl 5(x,) + VxS(x)} (26) 

in the ensemble members Xj, i = l,...,m. The equations are closed through the standard 
definition 



i=l 

for the mean and covariance matrix 

1 m 

P(s) = — I (Ms) ~ x(s)) (x,( S ) - x(s)) T . (28) 



i=i 



Since the covariance matrix P is symmetric, a straightforward calculation reveals that 



+ ( 29 ) 



d 
di 

along solutions of (126]) . More precisely, (126]) is equivalent to the gradient system 

^x, = -PV Xl V(X), (30) 
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in the ensemble matrix X(s) with potential 



m 



i=i ) 



V(X) = -<j5(x) + -^5(x l )^ (31) 



The actual decay of the potential V along solutions of (1301) depends crucially on the covariance 
matrix P. 

We note that Schur-product-based localizations can easily be applied to f )26|) to obtain, e.g., 



^-x, = -ip {V Xl S(x 4 ) + VxS(x)} , P = C loc o P, (32) 

or 

A Xi = -^(HPfR- 1 {Hx, + Hx - 2y} , HP = C loCjl o (HP), (33) 
ds 2 

in case of linear observation operators. These modified ensemble update formulations are easy 
to implement numerically. See Section H] for details. 

4 Numerical implementation aspects 

The various ODE formulations for the ensemble members Xj, % = 1, . . . , m, need to be solved 
over a unit time interval with initial conditions provided by the forecast values Xjj of the 
ensemble members. We apply the forward Euler method with step-size As = 1/4 (four time- 
steps) for our experiments. We found that As = 1 (single time-step) and As = 1/3 (three 
time-steps) lead to unstable simulations, while As = 1/2 (two time-steps) leads to occasional 
instabilities for larger values of the ensemble inflation factor 5 in ffT8l) . On the other hand, 
increasing the number of time-steps beyond four did not change the results significantly. We 
also expect that four time-steps will generally be sufficient in practical applications unless 
observations strongly contradict their forecast values and large gradient values are generated 
in (|26|) . The same consideration can apply to simulations with large inflation factors S in (ITS]) . 
As a safe guard, one can monitor the decay of the potential (13"T|) along numerically generated 
solutions and adjust the step-size As if necessary. 

Note that the continuous formulations do not require matrix inversions/factorizations except 
for the computation of R _1 . The computational cost ofjocalization can be reduced even further 
by using the following approximation. The matrix HP in f )33|) varies along solutions and an 
approximative formulation is obtained by replacing HP(s) by its value at s = for all s > 0. 
This leads to a linear ODE in the ensemble members x, with constant coefficient matrix, i.e., 



d 1 



-X,; 



HP(0)) T R _1 {Hxi + Hx - 2y} . (34) 



ds 2 

Note that HP an d HP H T are sparse matrices for compactly supported filter functions 



( IGaspari and Cohnl . Il999f ). Numerical implementations of f[3~4"j) should first update the incre- 
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Figure 1: The best RMS error for the Lorenz-96 model with an ensemble size of m = 10 and 
k = 20 observations taken in intervals of At b s = 0.05 over a total of 5000 assimilation cycles. 



ments Zj = Hxj — y with Euler's method, i.e., 




I = 0, 



(35) 



L = 1/As the number of integration steps, and then use the accumulated increments 

L-l 
1=0 



(36) 



to compute the final update of the ensemble members sq, i — 1, . . . , m. Overall matrix- vector- 
products will induce a computational complexity of O(km) in the ensemble size m and the 
number of observations k indepe ndent of the system size n. The same order of complexity 
applies to the serial algorithm of lHamill et al.l ( 1200 lh with the important difference that fl3T 
ca n be implemented as a si multaneous update over all observations. 



Bergemann et al.l ( 120091 ) proposed a re-orthogonalization technique for the ensemble devia- 



tion matrix X'. It should be noted that the re-orthogonalization is not uniquely defined. We 
implemented several variants of re-orthogonalization in combination with localization but did 
not find any significant improvements in the results. 



5 Numerical experiments 

We now report results from two test problems and implementations of (133]) and (134|) with local- 
ization. The results are compared to those from standard ensemble Kalman filter techniques. 
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5.1 Lorenz-96 model 



The standard implementation of the Lorenz-96 model ( ILorenzl . 1 19961 ; lLorenz and Emanuel! 



19981 ) has state vector x 
differential equations 



x, 



n 



40, and its time evolution is given by the 



x 



Xj + 



(37) 



for j — 1, ... , n. To close the equations, we define x_i = £39, xq = £40, and 241 = X\. 

The attractor of this standard implementation has a fractal dimension of about 27 and 
13 positive Lyapunov exponents. Localization will be necessary for ensembles with m < 13 
ensemble members. We use an ensemble size of m = 10 in our experiments. We observe every 
second grid point, i.e., k = 20, and the measurement error covariance is R = l2o- Measurement 
are taken in time intervals of At obs = 0.05. After a short spin-up period, a total of J = 5000 
analysis steps are performed in each experiment. The "true" trajectory x tmt h(t n ) is generated 
by integrating the Lorenz-96 model with the implicit midpoint rule and step-size At = 0.005, 
i.e., we assume that there is no model error. The observations are obtained according to 



y(^obs) — Hx tmth (t obs ) + r(t obs ) 



(38) 



where r(/; bs) are i.i.d. Gaussian random numbers with mean zero and covariance matrix R. 

We implement localization combined with standard ensemble inflation fo r the following 
five different ensemble Kalman filters: (i) EnKF with perturbed observations (IBurgers et al 



19981 ; iHoutekamer and Mitchell!. 120051). (ii) ensemble sq uare root filter (ESRF) with sequentia l 
treatment of observations (IWhitaker and Hamilll . |2002| ). (iii) DEnKF (ISakov and Okd . l2008al ). 
(iv) formulation (133]) . denoted CEnKF-I, (v) formulation (134)) . denoted CEnKF-II. We imple- 
ment CEnKF-I with As = 1/4 and As = 1/6, respectively, to demonstrate the impact of the 
discretization parameter on the results. 

For simplicity, localization is performed by multiplying each element of the matrices HP 
and HPH T , respectively, by a distance dependent factor pjy. This factor is defined by the 
compactly supported localization function (4.10) from lGaspari and Cohnl (119991 ) . distance r^i = 
min{|z — i'\,n — \i — where i and i' denote the indices of the associated observation/grid 
points Xi and Xy, respectively, and a fixed localization radius tq. The localization radius is 
varied between r = 2 and r = 30. The inflation factor 5 in (1181) is taken from the range 

[VTm,VTJE]. 

In Figure [TJ we display the RMS error 



rmse 



\ 



1 

— - } ||x(j • At obs ) - x truth (j • At obs ) 



(39) 



j'=i 



for an optimally chosen inflation factor 5 as a function of the localization radius r . Results are 
displayed for those localization radii tq only, which lead to at least one simulation with a RMS 
error of less than one. 

We conclude that EnKF yields the lowest filter skills while all other methods show an almost 
identical performance. 
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Table 1: Mean RMS error for localized CEnKF-I over 4000 time steps as a functions of the 
localization radius r and the inflation factor 5. For clarity, the value Inf is assigned if the RMS 
error exceeds the value 2.0 (no filter skill). 
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Table 2: Mean RMS error for localized CEnKF-II over 4000 time steps as a functions of the 
localization radius r and the inflation factor 5. 
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Table 3: Mean RMS error for localized DEnKF over 4000 time steps as a functions of the 
localization radius r and the inflation factor 5. 
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Figure 2: The best RMS error for the QG-model of ISakov and Okd (|2008aj) with an ensemble 
size of m = 25 and k = 300 observations taken in intervals of At Q b s = 4.0 over a total of 1000 
assimilation cycles. 



5.2 A quasi- geostrophic (QG) model 

We use the QG model of ISakov and "Okelfl2008af ). The QG model is a numerical approximation 
of the following 1.5-layer reduced gravity quasi-geostrophic model with double-gyre wind forcing 
and biharmonic friction: 



-ip x - eJ(ip, q) - AA 3 ip + 2irsm(2wy), 



(40) 



where q = Atp — Ft/), J(i/>, q) = i/j x q y — ip y q x , A = d 2 /dx 2 + d 2 /dy 2 . The coefficients are given 
by F = 1600, e = 10~ 5 , A = 2xl0~ 12 . The model domain is (x,y) G [0, 1] x [0, 1] with zero 
Dirichlet boundary condit ions. The mod e l is dis cretized over this domain using a 129 x 129 
grid. For more details see ISakov and Okd (j2008al ). 

We implement the deterministic ensemble Kalman filter (DEnKF) and our ensemble Kalman 
filters based on f )33|) . All experiments use m = 25 ensemble members. The dimension of the 
phase space is n = 16129. The dimension of the attractor and the number of positive Lyapunov 
exponents are c urrently not known. 



In line with ISakov and Okd (j2008al ). localization is performed by multiplying each element 
of the matrices HP and HPH T , respectively, by a factor Pi^vy = exp(— 0.5^ ^-,/^). Here 

we use the distance r^^f = — i'\ 2 + \j — j'\ 2 , where (i, j) and (i', f) denote the indices of 
the associated observation/grid points and a^y, respectively, and r is a fixed localization 
radius. 

We test the performance of the filters for different values of the ensemble inflation factor 5 
in ( ITS]) and the localization radius rr> For each pair (5, rn) of simulation parameters, we run 
a single simulation over 4000 time steps with step-size At = 1.25 and perform a total of 1000 
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assimilation cyc l es usin g 300 observations of ip with observation variance of 4.0 as described in 



Sakov and Okd (j2008al ) (after a spin- up period of 200 time steps and 50 assimilation cycles). 



All simulations are started from the same initial ensemble and use i dentical sets of observat ions. 

Since the DEnKF has been compared with EnKF and ESRF bv lSakov and Okd (j2008al ) and 
showed the best performance of all tested methods for this test problem, we only perform a 
comparison between the new formulations CEnKF-I (based on (13B1 with As = 1/4), CEnKF-II 
(based on ( I34p with As = 1/4) and DEnKF. The mean RMS errors of all three methods can 
be found in Tables [T] to In Figure [21 we display the RMS error for an optimally chosen 
inflation factor 5 as a function of the localization radius r$. Curves are based on the data 
presented in Tables [I] to [3j We conclude from Figure [2] that all three filters display a nearly 
identical performance for an optimal (5, ro) parameter choice and that DEnKF shows a slightly 
better p erformance for t he larg est localization radius r = 35. A similar observation was 
made bv lSakov and Okd (j2008af ) with regard to a comparison between DEnKF and ESRF. The 
differt results for r n = 35 could be due to the built-in overestimation of the analyzed ensemble 
covariance matrix (ISakov and Okd . l2008al ) . 

It should be noted that CEnKF-II is the least computational expensive of the three methods 
considered in this study. 



6 Conclusions and further extensions 

Schur-product-based localization of covariance matrices has become a popular and powerful tool 
to make ensemble Kalman filters perform well even under small ensemble sizes. In this note, 
we have proposed a new approach to implement Schur-product-based localization seamlessly 
within the framework of ensemble Kalman filters. Our approach is based on the formulation 
of the Kalman update step as differential equations in terms of its ensemble members. We 
have demonstrated for the Lorenz-96 model that the resulting methods outperform EnKF 
with perturbed observations and perform as well as standard implementations of ensemble 
Kalman filters such as ESRF with serial processing of observations and the recently proposed 
DEnKF. We also implemented a QG model and found that our methods perform nearly as 
well as DEnKF which is currently the best available method for this model problem. From 
a computational point of view, the formulation (134T) is particularly appealing since it leads to 
very efficient implementations without the need of matrix inversions (except when the error 
covariance matrix R is not diagonal) and only a single evaluation of the ensemble generated 
covariance matrix HP. 

We now outline an number of possible extensions of the formulation (13"01 . 

First we note that (|2T)1) can be used in connection with any cost functional 5(x) and, hence, 
provides a straightforward generalization of EnKF to nonlinear observation operators y = h(x), 
i.e., 

5(x) = i(h(x)-y) T R- 1 (h(x)-y). (41) 

Second, as for other localization techniques, the formulation (|32|) leads to updates in the 
ensemble deviations which lie outside the space S/ in general and, hence, may introduce 
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imbalance into the analyzed ensemble members Xj. It seems feasible to restore balance within 
the proposed framework by introducing additional cost functions S pS eudo( x ) the formula- 
tions f]3"2"]) or (1551) . respectively For example, we might require that the divergence of a velocity 
field v remains "small" by including a cost functional 



dV, 



*Spseudo 

where r > is an appropriate constant. Hence we would modify (125]) to 

5'pseudo(Xi)} + Vx {£(x) + Spseudofx)}] . 



(42) 



(43) 



Third, we have focused on deterministic ens emble Kalman filter formulations in this paper. 
However, EnKF with perturbed observations (IBurgers et all [1998) can also be put into the 
framework of continuous updates and leads naturally to the formulation 



—x, = -PH T R 1 {Hx, - y,} 



(44) 



where y« are now stochastically perturbed observations (IBurgers et all Il998l ). Alternatively, 
we may consider the stochastic differential equation 



dxi = -PH T R 1 {Hxids - yds + R 1/2 dwJ 



(45) 



in the ensemble members, where W j(s) G M fc denotes standard £;-dimensional Brownian motion. 
See, for example, I Gardiner! (120041 ) for an introduction to stochastic differential equations. 

Fourth, we have extensively discussed the continuous formulation of a single Kalman filter 
analysis step for a set of observations given at some time instance tj. We now come back to the 
complete ensemble Kalman filter formulation for sequences of observations at time instances 
tj, j = 1, . . . , M, and intermediate propagation of the ensemble under the dynamics (Q. The 
continuous formulation of the ensemble Kalman filter step allows for the following concise 
formulation in terms of a differential equation 



M 



/( Xi )-5>(*-^)PV Xi V;(X) 



(46) 



in each ensemble member, where S(-) denotes the standard Dirac delta function and Vj(X) is 
the potential (13T1) with 5(x) replaced by 



S,-(x) 



- (Hx 

2 v 



yfoOfR- 1 (Hx-yfe)). 



(47) 



One may view (146]) as the original ODE (Tjj) driven by a sequence of impluse like contributions 
due to observations. Numerically, it makes sense to regularize these impulses and to replace 
by 



M 



/(xO-^^-^PV^X), 



(48) 
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where 

ip(s) is the standard hat function (or a mollifier in the sense of iFriedrichs and e > is 

an appropriate parameter. Formulation (148]) can be solved numerically by any standard ODE 
solver. There is no longer a strict separation between ensemble propagation and filtering. Of 
course, the ODE (I48p becomes extremely stiff as s — > 0. A sensible choice is e ~ At, where At 
is the natural step-size for the ODE ([T]). The mollified formulation (148 p might be of particular 
interest in the context of the assimilation of no n-synopti c meas urements, e.g., measurements 
which are arriving semi-continuously in time (see 
approaches) . 



Evensenl (120061 ) for a discussion of alternative 
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